Nanoscale buckling deformation in layered copolymer materials 
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In layered materials, a common mode of deformation involves buckling of the layers under tensile 
deformation in the direction perpendicular to the layers. The instability mechanism, which operates 
in elastic materials from geological to nanometer scales, involves the elastic contrast between different 
layers. In a regular stacking of "hard" and "soft" layers, the tensile stress is first accommodated by a 
large deformation of the soft layers. The inhibited Poisson contraction results in a compressive stress 
in the direction transverse to the tensile deformation axis. The "hard" layers sustain this transverse 
compression until buckling takes place and results in an undulated structure. Using molecular 
simulations, we demonstrate this scenario for a material made of triblock copolymers. The buckling 
deformation is observed to take place at the nanoscale, at a wavelength that depends on strain rate. 
In contrast to what is commonly assumed, the wavelength of the undulation is not determined by 
defects in the microstructure. Rather, it results from kinetic effects, with a competition between 
the rate of strain and the growth rate of the instability. 

The mechanical response of multiblock copolymers depends sensitively upon their constituent homopolymers seg- 
ments, molecular architecture and chain topology. Triblock copolymers, in particular, have become an attractive 
material for their use as thermoplastic elastomers that could be integrated in several technical and manufactural 
fields ( copolymer styrene butadiene rubber is commercially exploited in footwear, in pressure sensitive adhesive 
(K-Resine), in paving and roofing compounds....). Depending on the amount of each phase the segregated block 
copolymers may present several morphologies e.g. spherical, cylindrical, lamellar... [1 j. The lamellar morphologies 
are particularly interesting as model systems, as well aligned specimens can be prepared by shearing. The one di- 
mensional aspect of the structure simplifies the analysis, and the mechanical response reveals equally the presence 
of both phases, while in other morphologies it tends to be dominated by the majority matrix phase. In these block 
copolymer systems the constituent blocks are generally chosen such that one of them is glassy and the other one is 
rubbery. A single copolymer chain can be shared between two different glassy lamellae, forming a rubber bridge that 
provides a strong coupling between phases. The resulting system combines the stiffness of the hard glassy phase and 
the ductility of the soft rubbery phase. 

When such a lamellar copolymer sample is submitted to a tensile strain perpendicular to the layers, the glassy layer 
eventually buckles into a "chevron" morphology. With increasing strain the normal to the lamellae tilts away from the 
stretching direction, whereas the lamellar spacing remains almost constant. This behavior was demonstrated experi- 
mentally in triblock copolymers by Small Angle Xray Scattering (SAXS) under deformation [2], and by micrographs 
of strongly deformed samples [3] (see figure [TJ. Similar effects are also observed in lamellar systems with alternating 
crystalline and glassy parts [1J[5], and in deformed semicrystalline polymers, which form locally lamellar structures 
of crystalline material separated by softer amorphous parts [6]. 

This buckling instability under strain, which is observed in many layered materials from smectic liquid crystals 
[7] to geological layers [9] , is frequently described in a qualitative way by a preference to shear compared to an 
extension in the direction normal to the layers, in order to preserve the lamellar spacing. A different cause for buckling 
is the existence of a Poisson effect. In a stacking of "hard" and "soft" layers, the soft phase accommodates most of 
the imposed deformation, and as the Poisson contraction is prevented by the coupling to the hard phase; it exerts 
a compressive stress in the transverse direction. The "hard" layers sustain this transverse compression until buckling 
takes place and results in an undulated structure. In general, elasticity predicts buckling to take place on the largest 
wavelength compatible with the boundary conditions imposed to the system [10] . We study this generic scenario by 
means of molecular dynamics simulations, for a material made of triblock copolymers in their lamellar phase. The 
contrast in elasticity is provided by a different glass transition temperature of the different blocks. We use the ability 
of molecular dynamics simulations to give information on the local values of stresses and strains to explore the causes 
of the instability in triblock copolymers with alternating glassy and rubbery layers, without introducing an a priori 
description of the mechanism. While the elastic origin of the instability is confirmed, our results demonstrate an 
unexpected dependence of the failure mode on strain rate. The wavelength of the undulations is determined by the 
deformation rate, and a characteristic size emerges even in the absence of preexisting defects in the microstructure. 
We discuss a general mechanism for the emergence of such a characteristic size, which results from a competition 
between the rate of strain and the growth rate of the buckling instability. 
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I. EVIDENCE OF THE POISSON EFFECT 

When a multilayered system is stretched perpendicularly to the layering direction, each component of the system 
will deform according to its own stiffness. Locally, the deformation is distributed between phases in a way that ensures 
the continuity of the stress. The resulting macroscale deformation is the sum of the local strain response of each phase, 
so that the stiffness will be dominated by the response of the soft phase. In a lamellar copolymer with alternating 
rubbery and glassy phases, the tensile strain will be mainly localized in the rubbery lamellae, which are essentially 
incompressible (Poisson ratio v ru bbery — 0.5). The contraction in the transverse direction that would result from the 
Poisson effect is however prevented by the strong coupling to the hard phase at the interface, so that the soft phase 
exerts a strong compressive stress in the transverse direction on the hard layers. The glassy phase becomes submitted 
to a tensile stress in the perpendicular direction and a compressive stress laterally. Under these conditions, and for a 
sufficiently large system, the buckling instability takes place to relax the lateral compressive stress. 

This scenario can be checked by monitoring the local stress in a sample strained perpendicularly to the lamellae. 
Figure [2] shows a map of the lateral stress ( ^ xx +^yy ) Q f a stretched copolymer sample at a true strain e zz = 0.04. The 
average pressure in the transverse direction is zero, as uniaxial tensile conditions are imposed globally. The local stress, 
however, is positive in the rubbery phase while it is negative in the glassy phase, indicating a local compression parallel 
to the glassy lamellae. An interesting additional feature that is apparent in figure [2] is the existence of noticeable 
fluctuations in the local stress values. Such fluctuations are expected in the response of amorphous materials, and 
arise from the local heterogeneity of elastic properties in such materials, a property already well documented for glassy 
polymers [TT] . We will see below that the response of the system can become macroscopically non uniform due to the 
development of an elastic instability, and it is therefore quite natural to speculate that the initial stress fluctuations 
serve to nucleate this instability. A related study of cavitation in glassy homopolymers showed a similar correlation 
between elastic fluctuations and nucleation of cavities [T2] . 

II. PLASTIC DEFORMATION THROUGH BUCKLING AND CAVITATION 

Figure [3] illustrates the deformation process that takes place under tensile deformation in a sufficiently large sample, 
and the corresponding stress-strain curve is shown in figure [4] In this first set of results, the strain rate is e zz = 
7.3 x 10 -5 in reduced units (see Materials and Methods section) . The pressure is fixed to zero in the X and Y 
directions. 

The stress strain curve displays three main regimes: 

• Elastic regime: at small strain, the stress grows linearly with the imposed deformation. This regime is 
limited to a very small deformation amount (less than 2%). The Young modulus, which can be fitted from 
this curve, results from a combination of the elastic behavior of each phase. The density maps show very little 
change in this regime. 

• Buckling : Beyond the elastic regime, a progressive softening is observed. This slight deviation from the elastic 
linear behavior is commonly interpreted by the change of molecular conformation in the rubbery phase. The 
buckling of the glassy phase starts at e zz — 0.06 and can be detected from a rapid change in the apparent 
Poisson ratio of the material. From the onset of the buckling instability, the sample starts deforming in an 
"accordion" like manner, with a higher Poisson ratio than a sample with parallel layers. 

• Cavitation: After buckling, a strong stress drop occurs at e zz = 0.08. This drop can be correlated with the 
nucleation of cavities in the rubbery phase, as illustrated by the fourth density map in figure [3] The low density 
spots in figure [3] correspond to the cavitation in the rubbery layers. Indeed, due to the buckling, the local 
deformation of the rubbery phase is not homogenous. The sample progressively adopts a chevron morphology, 
with different deformation states: at large strains the cavitation is confined to a localized region in space, 
developing into a hinge. At the hinges of the chevron the deformation is essentially tensile, while the tilted part 
undergoes a simple shear deformation. The latter is caused by the rotation and sliding of the hard lamella. In 
contrast, at the hinges of the chevron the deformation is essentially triaxial, and favors nucleation of cavities. As 
a result the cavities in the rubber that initially appear randomly tend to heal in the sheared zones and nucleate 
preferentially where triaxial stress persists, as illustrated by the sequence of snapshots in figure [3] ). 
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FIG. 1: Top panel: Transmission electron micrograph shows the buckling observed in a real copolymer sample (Styrene- 
Butadiene-styrene) . The figure is extracted from reference [13] . Bottom panel: Lamellae buckling observed by molecular 
dynamics simulation, snapshot is taken at e zz — 0.55. This image contains six periodic copies of the simulated system, which 
is enclosed in the white frame. The scale on the right is obtained by assuming a monomer size a = 0.5nra. 

III. ELASTIC DESCRIPTION OF THE BUCKLING INSTABILITY 

As mentioned in the introduction, two different theoretical descriptions of buckling of lamellae under stress are 
available in the literature. One approach is based on a free energy functional of the order parameter that describes 
the lamellar order, and explains the instability by the fact that a strained state will try to maintain the lamellar 
distance that minimizes this free energy [14]. This approach would be appropriate for copolymers in which both 
phases are at equilibrium, so that the order parameter can respond to the deformation. Here one of the phases 
is glassy, and the explanation for the buckling instability must be searched in a different direction, involving the 
minimization of the total elastic energy. This description (which for some aspects goes back to early works of Biot 
[H[T5j[T6]) has been detailed in a seminal paper by Read et al [ID] . Taking for simplicity a two dimensional geometry, 
the elastic energy of the sample under small strain can be written as: 

Umacro = ^(Ciien + 2C 12 ene 2 2 + C 22 el 2 + Ge 12) (!) 

where Cij are the components of a symmetric 2x2 stiffness matrix of the entire system, and e^- is the macroscopic 
deformation. In this equation the material is considered as homogeneous and anisotropic. The indexes 1,2 refer to 
the directions normal and parallel to the layers, and the coefficients Cij and G (shear modulus parallel to the layers) 
can be obtained from mechanical tests in different directions, or deduced from the elastic properties of each lamella. 

The bending energy results from the variation of the lamellar rotation angle (the angle between the layer and the 
horizontal axis X ). The associated energy density reads 



Uhend = \k(V x 9) 2 



(2) 
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FIG. 2: Cartography of the lateral stress xx 2 — — . Snapshots were taken at a deformation e zz — 0.04. The loca l str ess is 
computed within a box of size 5cr x 5cr x L y — 102cr, where a is the size of a monomer (typically a — 0.5nra, see section VI ). The 



stress in these figures was probed in a relatively small sample (L x = 34.5cr, L y = 102. 3a and L z = 34.5cr). in which buckling 
does not take place, however the analysis of the stress profile clearly shows the importance of the Poisson effect, with a large 
contrast in the transverse stress between the compressed glassy layers and the rubbery layers under tension . The influence of 
the lateral dimension of the sample is discussed below. 
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FIG. 3: Local density cartographies at several strains. The high (low) density phase corresponds to the glassy (rubbery) 
lamellae. The phase buckling starts before cavitation. As the deformation progresses, cavities will nucleate randomly in the 
rubbery phase. Cavities that are located in the tilted part of the chevron disappear rapidly; however, only cavities that are 
located in the hinges will survive to a high strain. 
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True strain e zz 

FIG. 4: Top panel: stress strain curve of the sample shown in figure [3] under uniaxial strain conditions. The black points denote 
the selected configurations shown in figure [3] The stress release is associated with a combination of two phenomena: buckling, 
and nucleation and growth of cavities in the rubbery phases. Bottom panel: evolution of the lateral strain 0.5(e xx + e zz ) under 
uniaxial tensile conditions. The local slope defines the Poisson ratio, and buckling is correlated to a strong change in this slope, 
as the sample starts deforming in an accordion-like manner. 



where K is the bending modulus of the sample. Due to the serial coupling between phases the bending modulus will 
be dominated by the contribution of the hard phase. The bending modulus can then be estimated from simple beam 

bending theory as K est = f h ( f hd 2\ where (j>h is the volume fraction of the hard phase, Eh is the Young modulus of 

1^1 v h ) 

the hard phase, Uh is its Poisson ratio and d is the lamellar spacing. 

The total energy density results from the addition of the bulk elastic energy and the bending energy, U2D = 
Ubend + Umacro. It can be expressed in terms of the global deformation e xx , e xz and e zz and of a local, non affine 
displacement field u(x, z) that quantifies the difference between the displacement imposed at the global scale and the 
one observed locally. 

The final expression for the free energy, written as an expansion in powers of u(x, z) [10 reads 

2(U 2d ) = Cne 2 xx + 2C 13 e xx e zz + C 33 e 2 zz 
+Cn{(V x u x ) 2 ) + 2C ls ((y z u x )(y x u z )} + C 33 ((\7 z u z ) 2 ) 

+ a+tlP \° ~ e -( C 33 - C13 - 2G) - e 2 zz (C 33 - G) (3) 
+e xx (C n (l + e xx ) - C l3 )\ + K{(V 2 x u z ) 2 ) 

+2 ((V^)^)) [G _ ezz{Cs3 _ Q) _ exxCis] + 0(m4) 

A linear stability analysis of this energy is performed by introducing a sinusoidal perturbation of the form observed 
in our simulations, 

u z (x, z) = Uq sm(kx) ; u x (x,z) = (4) 
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with k = 2nw/L a wave vector compatible with the boundary conditions, 

2(U 2 d) = C u e 2 xx + 2C 13 e xx e zz + C 33 e 2 zz 

(5) 

+ \{fi(e xx , e zz )k 2 + Kk 4 }U 2 + O(C/ 4 ) 

where 

fi D (e xx ,e zz ) = )2 [G - e zz (C 33 - C 13 — 2G) 

(6) 

-e 2 zz (C 33 -G) + e ra (Cn(l + e xx ) - C 13 )} 

The buckling instability occurs upon increasing strain when the coefficient of U§ in equation [5] becomes negative, 
meaning that the global gain in elastic energy overwhelms the bending energy penalty. In 3D case the same analysis 
process can be performed, this leads to a function fi D (e xx1 e yyi e zz ) instead fl D (e xxi e zz ) where: 



fl &zz-> Cyy) — fi i^xx^zz) n j r&xx ^2 ClZ^yy (7) 

To close the system, one assumes that before the buckling begins (i.e. in the elastic regime) e xx and e yy can be 
substituted by ve zz where v is a global Poisson ratio. Under this assumption, f\ D (e xx , e yy , e zz ) becomes a function 
fi D ( e zz) of e zz only. The buckling strain e\ uck can be estimated by solving this equation for a fixed wavevector 
k n = 2rmjL\ 

fi D (eLck)kl + Kk*=0 (8) 

For a given wavevector k n , buckling will become possible above a certain strain e* such that /| D ( e n) = —Kk\. As 
| /I" 1 is an increasing function of the strain, the wavevector corresponding to the largest wavelength, i.e. the size of 
the box, will become unstable at the smallest strain, according to this analysis. This also implies that for a smaller 
box size, a larger strain would be needed to observe buckling; as noted above, in such small systems cavitation tends 
to take place before the critical strain for buckling is reached, and the elastic analysis becomes irrelevant above the 
cavitation threshold. 

In the following, numerical comparison between simulations and this theory will be made by using for the elastic 
constants Cij values determined from simple linear deformations of a small sample that does not exhibit the buckling 
instability. These values are, for the interaction parameters and temperature mentioned below, C\\ = 24.17, C33 = 
7.61, C23 = C13 = 6.5, and G = 0.07 (in the same reduced units, the Young modulus of a glassy polymer is of order 
50, see section "Materials and methods"). The Poisson ratio is v = 0.178. 



IV. SIZE AND STRAIN RATE DEPENDENCE, INTERPRETATION 

The discussion in the previous section shows that the boundary conditions and geometry of the sample have an 
important influence on the buckling instability. For example, the instability takes place for a different mode for a 
system with free or with periodic boundary conditions, the first one undergoing a "half wave" instability which is 
prohibited in the second case. Also, the buckling of small samples is impossible as the bending energy of the glassy 
phase is very high compared to the deformation energy of the bulk. Therefore, a critical size of sample can be defined 
as L*. L*, the minimal length from which the sample will be able to buckle under tensile strain before failure due to 
cavitation and crazing in the rubber phase becomes dominant. 

According to elastic theory, the instability will take place at smaller and smaller strains for bigger and bigger 
samples, and always at the largest possible wavelength allowed by the boundary conditions. This, however, contradicts 
a number of experimental observations in which a rather well defined wavelength of the chevron structure is observed. 
This discrepancy is commonly assigned to preexisting defects in the microstructure [3]. Simulation provides an ideal 
benchmark of this hypothesis, as it allows one to study an ideal microstructure. We have therefore studied samples 
of various sizes, between 100 and 800a in the Y direction. Figure [5] compares the response of all tested samples, at a 
strain rate 7.3 x 10 -5 . In terms of stress-strain relation all samples have roughly the same mechanical behavior up to 
the yield point. 

However, it is clear that the minimal length for observing buckling before cavitation is at least 200 monomer sizes 
cr, implying that such observations can be made only on very large (by simulation standards) systems. For samples 
larger than about 400<r, the onset of buckling occurs always at the same strain (eb uc k = 0.06), in contradiction with 
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FIG. 5: Snapshots of several samples with different sizes. All configurations are taken at a true strain e zz — 0.08. The samples 
are obtained by replication of a system with size Lo = 32.8cr. Buckling is possible only for the longer samples with L y > 393.6). 
The buckling wave length seems to be size independent. 
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FIG. 6: Strain at buckling. Comparison for two different strain rates 7.3 x 10 -5 (blue bars) and 1.4 x 10 -5 (red bars) for 
different system sizes, with the prediction of elastic theory (yellow bars). Note that the smallest system does not exhibit 
buckling at high strain rate. For large samples, the buckling at high and low strain rates are observed in different modes, see 
figures [7] and [8] The comparison with the prediction of elastic theory shows that the latter is quite accurate when the buckling 
takes place in the largest wavelength mode, although it underestimates somewhat the buckling strain due to kinetic effects. 

the expectation from the elastic description of the previous section. Another surprising observation, illustrated in 
figure |5j is that the wavelength of the instability does not appear to increase with the size of the system. 

Figure [6] summarizes the buckling properties obtained at strain rate of 7.3 x 10 -5 . For each sample, the theoretical 
value of the buckling strain was calculated from equation [8j assuming an instability wavevector k = 2ir/L except 
for the largest sample. The predicted values of e b ^ z ck is always less than the measured one. This difference will be 
interpreted below as a direct consequence of kinetic factors that are not taken into account in the elastic calculation. 
Note, finally, that the yield strain and stress are roughly the same for all samples, and are correlated with the 
cavitation in the rubber phase. 

In order to understand the role of strain rate, we have submitted the same samples to a uniaxial tensile strain test 
driven at a lower strain rate (e yy = 1.4 x 10 -5 ). We find that (i) the change in the Young modulus is negligibly small 
for all samples, (ii) the yield stress and strain decrease as the strain rate decreases and finally (iii) the stress softening 
exhibits a smooth transition (from yield to the drawing regime) at low strain rate compared to a large drop at high 
strain rate. Depending on sample size, the yield stress and strain are more or less affected. For the smallest sample, 
the decrease of the yield stress and strain is small compared to other samples. In general, the decrease of the yield 
threshold is strongly correlated with the change of the plastic mode from cavitation to buckling. Both cavitation 
and buckling result in a yield behavior, however the yielding associated with buckling is much more progressive and 
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FIG. 7: Top panel: snapshots of samples with lateral size L y — 492a and L y — 787a driven at two different strain rates 
(e'yy = 1.4 x 10 -5 and e yy = 7.3 x 10 -5 ) . At low strain rate , the buckling wave length is equal to the sample length which 
is not the case at high strain rate. Note the imperfect aspect of the sinusoidal pattern in the first snapshot, indicative of a 
competition between modes with different wavelength. Note also that for this high driving rate, cavitation and buckling take 
place simultaneously. Bottom panel: stress strain curves that correspond to the two different driving rates of the largest sample. 
The high strain rate curve is identical to the one shown in figure [4] 

smooth than the one associated with cavitation, as illustrated in figure [7] The most striking observation, however, is 
the fact that the wavelength of the buckling instability changes with strain rate. In the largest sample, it is equal to 
sample size, as expected from elasticity theory (see figure [7]). 

We propose an interpretation of this unexpected change of wavelength with strain rate based on a simple model 
of the competition between unstable modes. It has been argued above (section " elastic description of the buckling 
instability") that the onset of buckling at a certain wavevector k is determined by the function F(e ZZl k) = fii^zz)^ + 
Kkf where is obtained from equation [6] or [7] Above a critical strain e c , /2(e) becomes negative, and the wavevectors 
in the range < k < \J—f2K become linearly unstable. We propose to account for the time evolution of the amplitude 
Uk of a given mode trough the simple linear equation 

dtU k = -AF(e zz (t),k)U k (9) 

where A is a phenomenological kinetic coefficient which is assumed here, for simplicity, to be wavector and strain 
rate independent, at least in th erange of strain rates investigated in our simulations. Beyond e c , the growth rate AF 
has a maximum at a finite wavevector. The key in understanding the effect of strain rate is in the time dependence 
of F, which is encoded in the time dependence of e zz and renders the evolution nontrivial. For small strain rates, 
the instability of the smallest wavevector k = 2tt/L has ample time to develop before the strain increases and makes 
a second mode at k = 4ir/L unstable. At higher strain rates, the growth rate evolves on time scales comparable to 
the instability itself, and a competition between modes arises. By solving numerically equation [9] with a value of 
A adjusted to reproduce the observed growth of the instability in one of the configurations, we found that, for the 
largest sample size considered in our simulations, the mode that corresponds to k = 4tt/L indeed overwhelms the long 
wavelength mode k = 2ir / L when the system is driven at the largest strain rate, before the latter mode has time to 
develop significantly. Therefore the second mode, rather than the first, is observed. While the very high strain rates 
used in simulation (see section "Material and method" ) result here in a selection of very short wavelengths of the order 
of nanometers, the effect is expected to be very generic, and could result in the submicron chevron structures observed 
in experiments. These results are summarized in figure [§J which displays the competing mechanisms depending on 
sample size and strain rate. 

V. CONCLUSION 

The buckling instability of layered materials has been explored here for the first time by atomistic simulations, in the 
case on a lamellar copolymer with glassy and rubbery lamellae. This approach allows one to identify unambiguously 
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FIG. 8: Qualitative summary of our observations, indicating the various modes of deformation and yield depending on strain 
rate and sample size. The points indicate the simulated configurations, while the solid lines between regions are shown as 
guides to the eye. 



the mechanism that drives the instability and failure of such materials, as being caused by the elastic contrast between 
the two phases and a Poisson ratio effect. Still, a pure elastic theory is insufficient to account for the dependence 
of the observed instability on strain rate, which is also reported here for the first time. A simple description of this 
dependence was provided by a model that involves the competition between the growth rate of linearly unstable modes 
with the rate of deformation. This rate dependence of the deformation mode of the microstructure, which is still to 
be investigated in controlled experiments, indicates the possibility of a rich and still unexplored phenomenology in 
the mechanical behavior of layered materials. The generality of the mechanism for buckling points to the importance 
of this instability as a new plasticity mechanism in nanostructured polymer materials, which is expected to operate 
in semicrystalline materials as well as in block copolymers, and to complement more traditional mechanisms such as 
cavitation and crazing. 



VI. MATERIALS AND METHODS 



Our molecular dynamics simulations involve a standard coarse grained model polymer chains. Each chain is repre- 
sented by a sequence of 50-100-50 beads of different chemical nature (A and£>) interacting via covalent FENE bonds 
and Lennard Jones interparticle potentials. The parameters of the LJ potentials are chosen to ensure that phase 
separation into a lamellar morphology takes place, and that the A layers are in a glassy state while the B layers 
remain rubbery at the tensile test temperature (T test ~ 300K). The method used for generating samples has been 
detailed in reference [17]. In terms of Lennard- Jones interaction parameters, we take caa = 1, £bb = 0.3, cab = 0.4. 
The resulting glass transition temperatures are c± 420K for the A layers and ~ 200K for the B layers. The 
diameter and mass, a and m, are identical for all species, and serve as length and mass unit, respectively. Using an 
energy scale of lOOOK/ks and a length scale of O.bnm which are typical in the coarse grained descriptions of standard 
polymers, the corresponding stress unit is of order lOOMPa, and the Young modulus of the glassy polymer is of the 
order of l-10GPa. 

An important remark here is that the time scale that results from these choices of units, if parameters appropriate 
for typical polymers are used, lies in the picosecond time range. Therefore the strain rates achieved in simulations 
are of the order of 10 T s -1 in real units, extremely high compared to typical experimental rates. As is often the case 
in simulation studies involving glassy materials, the behavior observed in simulation studies must be understood as 
being qualitatively, rather than quantitatively, representative of the experimental reality. 
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